Aging in short-ranged attractive colloids: A numerical study 
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We study the aging dynamics in a model for dense simple liquids, in which particles interact 
through a hard-core repulsion complemented by a short- ranged attractive potential, of the kind found 
in colloidal suspensions. In this system, at large packing fractions, kinetically arrested disordered 
states can be created both on cooling (attractive glass) and on heating (repulsive glass). The 
possibility of having two distinct glasses, at the same packing fraction, with two different dynamics 
offers the unique possibility of comparing — within the same model — the differences in aging 
dynamics. We find that, while the aging dynamics of the repulsive glass is similar to the one 
observed in atomic and molecular systems, the aging dynamics of the attractive glass shows novel 
unexpected features. 

PACS numbers; 64.70.Pf, 82.70.Dd 
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I. INTRODUCTION 

In the last years, the physics of systems characterized 
by short-ranged attractive inter-particle interactions — 
short in comparison to the particle size — has been the 
focus of several investigations. A good example of this 
type of systems are colloidal solutions, since the intensity 
and the range of the potential can be finely controlled by 
changing the solvent properties or the chemistry of the 
dispersed particles. The short range of the attractive 
interaction produces thermodynamic features which are 
not observed in atomic or molecular system (where the 
interaction range is always long ranged) [E S M ■ Re- 
cently, theoretical, numerical and experimental studies 
have focused on the dynamical properties of these sys- 
tems. One of the most astonishing discovery is that, at 
high density, the metastable liquid is characterized by a 
non monotonic temperature dependence of the charac- 
teristic structural times: the dynamics slows down not 
only upon cooling (as commonly observed in molecular 
systems) , but also upon heating. The slowing down upon 
heating can be so intense that a novel mechanism of ar- 
rest takes place high T 0. In the high density part of 
the phase diagram, for sufficiently short-ranged poten- 
tials, a re-entrant liquid-glass line is observed and a liquid 
phase emerges between two glasses. In other words, mov- 
ing along a constant density path, it is possible to pass 
from a glass phase to a liquid phase and then again in a 
new glass phase just by progressively lowering the tem- 
perature. The two glasses are named "attractive glass" 
and "repulsive glass" due to the different mechanisms, re- 
spectively attraction and excluded volume, which control 
the dynamical slowing down ^. In the re-entrant liq- 
uid region, between the attractive and repulsive glasses, 
dynamics is very unusual. Correlation functions show a 
subtle logarithmic decay, while the particle mean squared 
displacement follows a power-law in time. 

The complex dynamics in short-ranged attractive sys- 
tems has been interpreted within Mode Coupling Theory 



(MCT) Theoretical predictions do not depend on the 
detailed shape of the attractive potential and have been 
confirmed for r nany different modelization of the attrac- 
tive tail. H H Efl m El. A set of simulations [H El 
E1E3 and experiments [H El ESSISl El 

have con- 
firmed the existence of such complex dynamics. 
In numerical work based on the Square Well System 
(SWS), one of the simplest models studied, it has been 
shown that a diffusivity maximum in temperature ap- 
pears when the width of the square well is short enough 
compared to the hard core diameter E3 • The calculation 
of the iso-diffusivity lines in equilibrium E3 confirms the 
reentrant behavior and the novel dynamics which take 
place in the reentrant region. Along an iso-diffusivity 
path, the decay of the density-density correlation func- 
tions, which at high temperature can be well described by 
a stretched exponential, becomes more and more uncom- 
mon upon cooling, showing at low T signs of logarithmic 
decay [H El m. 

Previous numerical work has focused on equilibrium 
properties and on the dynamics near the predicted glass- 
glass transition 0, 0, 0, . Here we report a study 
of the non-equilibrium dynamics. Starting from an equi- 
librium state in the re-entrant fiuid region, the system is 
quenched both at low and high temperature. We are then 
able to explore — within the same model — the aging 
dynamics for the attractive and repulsive glass, to find 
out if the differences which characterize the dynamics in 
the glassy state show up also in the out-of-equilibrium 
evolution. We discover that, while the aging dynamics 
of the repulsive glass is similar to the one observed in 
atomic and molecular liauids[25ll2fi||. the aging dynamics 
of the attractive glass shows novel unexpected features. 
We also report comparison with a recent experimental 
study of aging in short ranged colloidal systems 
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II. SIMULATION DETAILS 

We investigate a system that has been extensively 
studied earlier in the equilibrium regime, a binary square 
well mixture [rsf. The binary system is a 50%-50% mix- 
ture oiN = 700 particles. The two species (labeled A and 
B) are characterized by a diameter ratio S — gaI'^b = 
1.2. Masses are chosen to be equal and unitary. The at- 
traction is modeled by a square well interaction defined 
according to: 

{oo rij < (Tij 
-uo Oij < nj < (Tij + Aij (1) 
r,y > o-y + A,y 

where cry = (ctj + crj)/2, i,j — A, B and has been 
chosen in such a way that e,;j = Aij/{dij + Aij) = 0.03. 
We choose fcs = 1 and the depth of the potential uq — 1. 
Hence T — 1 corresponds to a thermal energy equal to 
the attractive well depth. The diameter of the small 
specie is chosen as unity of length, i.e. as = 1- Den- 
sity is parametrized in terms of packing fraction (j) = 
{pAd\ + PBd%) ■ 7r/6, where pi = Ni/L^, L being the 
box size and iVj the number of particles for each species. 
Time is measured in units of as ■ (m/ugY^^ . The unit has 
been chosen consistently with Ref. [T^, where the same 
system has been studied in equilibrium. A standard MD 
algorithm has been implemented for particles interacting 
with SWS potentials Between collisions, particles 

move along straight lines with constant velocities. When 
the distance between the particles becomes equal to the 
distance where V{r) has a discontinuity, the velocities 
of the interacting particles instantaneously change. The 
algorithm calculates the shortest collision time in the sys- 
tem and propagate the trajectory from one collision to 
the next one. Calculations of the next collision time are 
optimized by dividing the system in small subsystems, so 
that collision times are computed only between particles 
in the neighboring subsystems. 

With the present choice of parameters, the formation of 
crystallites during the simulation run is not observed. 
As compared to the monodisperse case, the slight asym- 
metry in the diameters allows us to follow the devel- 
opment of the slow dynamics over several orders of 
magnitude We choose to work at the packing frac- 
tion (p — 0.608, being this density in the middle of the 
re-entrant region of the phase diagram, as shown in pre- 
vious calculations . At this density we prepared 60 in- 
dependent configurations, equilibrated at the initial tem- 
perature Ti = 0.6. Subsequently we quenched these inde- 
pendent configurations to the desired final temperatures, 
and then we followed the evolution in time at constant 
temperature. The characteristic time of the thermostat 
as been chosen to be much smaller than the structural re- 
laxation time and such that the system may equilibrate 
within one time unit. Each of the 60 independent con- 
figurations has been quenched to Tf = 1.2, i.e. in the re- 
pulsive glass, and to T/ = 0.3, i.e. in the attractive glass 
and run up to = 2 • 10^. The numerical protocol is 



illustrated in Fig. Iso-diffusivity curves from Ref. 
and the extrapolated ideal glass line from Ref. are 
plotted together with the quench path selected in this 
work. Since the system, after the quench, is out of equi- 
librium, all averages — represented in the text by the 
symbol (•) — have been performed over the ensemble of 
the 60 initial independent configurations. 



III. RESULTS AND DISCUSSION 

Numerical investigations of aging in glassy materials 
are usually accomplished by rapidly quenching the sys- 
tem from an equilibrium state to a new state, where the 
characteristic time scale of the structural relaxation is 
much longer than the observation time. When the time 
scale of the final state exceeds the characteristic time 
of the simulation, equilibrium cannot be reached within 
the simulation time. Out of equilibrium time transla- 
tion invariance does not hold and the state of the system 
does not only depend from the external parameters, but 
also from the time elapsed from the quench, the so-called 
waiting time tw ■ Indeed in the aging regime other notable 
effects e merge, i. e. violation of the fluctuation dissipation 
theorem |33, |^ Is^ l3^ and strong dependence on 
the story of the system |25l |. Thus, in out-of-equilibrium 
conditions correlation functions will depend both on the 
time of observation t and on the waiting time tu,, i.e. 
CA{tw,t) = {A* {tw)A{tiu Similarly, the one time 

properties (e.g. the static structure factor or the energy), 
will depend on the waiting time t^,. 
The first quantity that we investigate is the average po- 
tential energy of the system U. At = 0+, i.e. instanta- 
neously after the quenches to respectively Tf = 1.2 and 
Tf = 0.3, the two systems possess the same potential 
energy. Then, with tj^ , they start evolving toward differ- 
ent energy values. The situation is represented in Fig. |21 
While the system quenched to Tf = 1.2, i.e. in the re- 
pulsive regime, shows an increase in the potential energy 
(i.e. a progressive breaking of the interparticle bonding), 
the one quenched to Tf = 0.3 shows the opposite behav- 
ior. In both cases, the time dependence of the potential 
energy shows a fast short-time (t^ < 1) evolution (which 
reflects the equilibration of the fast degrees of freedom) 
followed by a much slower evolution, proper of the aging 
dynamics. The aging dynamics has detectable different 
trends at high and low temperatures. While at high tem- 
peratures a clear log{tw) dependence of U is observed, at 
low T the time dependence in log scale is not linear and 
can be well fitted hy U — Ueq t~°-, with a ~ 0.14. A 
similar form was successfully adopted by Parisi for a bi- 
nary system of soft spheres in the aging regime [30l| . 
The time evolution of the aging structure can be followed 
by numerical calculation of the partial static structure 
factors, 5*0/3(9, tto) defined as 

Sal3{q,tw) = {g*a{<i,tw)gp{<l, tw) / y/naTlf}) (2) 
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where the partial density variables are defined as Qa {q) = 
exp[iq • r(i)^"'']/\/]V, and rii — Ni/N are the partial 
concentrations. In equilibrium, time translation invari- 
ance implies Sai3(q,t) = Saf3{q,0), i.e. the static struc- 
ture factor is a time independent function. 
The inset of Fig. |31 shows SAA{q,tw) for the quench to 
Tf = 1.2, for waiting times spanning 4 decades. Similarly 
to binary Lennard- Jones systems j2a| , the variation of the 
static structure factor is particularly small, for all a, f3 
components. With increasing a very weak increase 
at all peaks of the structure factors is observed. This 
is shown in the main figure. Similarly for the Tf — 0.3 
case, the Sap{q,tw) changes are very small. Interestingly 
enough, for the first peak (and only there) the trend is re- 
verted with respect to the previous case, i.e. the intensity 
of SAAiq,tw) around the peak decreases on aging. This 
observation agrees with equilibrium studies. Indeed, for 
the square well system, an increase in temperature pro- 
duces an increase in the first peak toward the hard sphere 
limit whereas lowering the temperature only promotes 
the structure factor oscillations at large wave-vector^S^. 
It is perhaps interesting to recall that within MCT, the 
increase of Sa/s at the first peak is responsible for the 
formation of the repulsive glass, whereas the long lasting 
oscillations in the structure factor generates the kinetic 
arrest in the attractive glass |T3|. 

Next we study the intermediate scattering function, 
i.e. dynamic density-density correlators. For a binary 
mixture this quantity is defined by: 

<Pap{q, t, t') = (ea(q, t)*gp{q, t'))/Sap{q, t') (3) 

where the density variables ea(q, i) have been defined 
before. In equilibrium 4>ap{,q,t,t') will depend only on 
the difference t — t' . When the system ages the density- 
density correlators is a function of and t. In this pa- 
per, we show (papiq, t, t') as a function of and the time 
difference {t — tw), i.e. 4>apiq,tw,t — tw). The density 
correlation is normalized by Sapiq, t^). The upper panel 
of Fig. [SI shows the dependence of — t^), 

with quB = 8.72, which lies between the first and the 
second peak in the static structure factor. On increas- 
ing tw, the decay of the correlation functions becomes 
slower and slower. The time needed to lose memory of 
the starting configuration grows. This is a typical feature 
of aging systems and it has been encountered both in ex- 
periments |27| and simulations pd[ . The slowing down 
of the dynamics on increasing resembles the way the 
structural arrest emerges in supercooled liquids on ap- 
proaching the glass transition. Indeed, in this case, upon 
cooling two characteristic time scales emerge, a fast de- 
cay around a plateau value and a slow decay to zero. The 
time duration of the plateau becomes longer, the closer 
the system is to the glass transition. On the other hand, 
in the aging regime the control parameter is the waiting 
time tu,, which plays a role analogous to that played by 
T in equilibrium, i.e. the longer 1^ the longer is the time 
requested to lose memory. In the same way as temper- 
ature measures the distance from the glass transition in 



equilibrium liquids, measures the progressive thermal- 
ization of the structural degrees of freedom. In summary, 
with increasing t„ the system develops a clear plateau, 
whose height is consistent with the value expected for the 
repulsive glass. 

Similar results are shown in the lower panel of Fig. for 
qaB = 6.37, which is close to the maximum of the static 
structure factor. Two main differences emerge between 
the correlators at the the two wavevectors: (i) the relax- 
ation time is longer for qaB — 6.37, (ii) the plateau is 
higher. Both this issues follow, as expected, the trend 
of the equilibrium results, where at the maximum of the 
structure factor a maximum of the relaxation time and 
of the non-ergodicity parameter is encountered j35j | . 
The upper panel of Fig. shows the density-density 
correlation functions for the quench at Tf = 0.3 for 
qcjB = 8.72. Similarly to what we found in the previ- 
ous case, the time of decay of the correlators grows with 
the waiting time t^. However the shape of the functions 
is drastically different, confirming that even the aging 
dynamics reflects the differences noted in the previous 
equilibrium studies. Two notable differences with the 
Tf = 1.2 case are: (i) the strength of the Q;-relaxation(i.e. 
the height of the plateau) in aging is much more intense, 
being now close to one ; (ii) the plateau does not sig- 
nificantly stretch in time (as shown in the inset) beyond 
t — tw ~ 1. Similar trends arc found for a g- vector closer 
to the first peak of the structure factor, as shown in the 
lower panel of Fig. El 

A recent study of the glass-glass transition has sug- 
gested the possibility that the attractive glas s is signif- 
icantly destabilized by hopping processes [2j, which in 
the present context can be associated with the breaking 
of the interparticle bonds. The data shown in the inset of 
FigEl support such interpretation, since no clear plateau 
can be observed for > 1. 

A very interesting feature of the correlation functions 
in the aging regime is the property of scaling with wait- 
ing time. Theoretical studies suggest that correlators at 
different waiting time can be rescaled by the equation: 

Kp{q,tw,t-tw) = (t)%{q,t) + (i)ap(q, ^^^^^^y^) (4) 

where (f'afsil^^) accounts for the short time (s^. The 
h{t) — t case is the so-called simple aging case. For 
a binary mixture of Lennard-Jones spheres the system 
follows gl) with h{t) = t" and a « 0.88 
We tried to apply this scaling relation to our system. 
We focus on the low temperature quench, since the high 
temperature one shows a phenomenology, which is similar 
to the Lennard-Jones case. In the case of aging of Tf — 
0.3 no clear scaling of the correlation functions can be 
performed. In Fig [71i we show the crude tentative of 
data scaling. Correlators have been shifted in time to 
superimpose the correlators in the region where 0.7 < 
4>q < 0.9. Fig. [7}p shows the same data, this time shifted 
to maximize superposition for 0.35 < (pq < 0.5. In Fig. |S1 
the scaling coefficients tr used to obtain the scaling in 
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Fig. 13 are shown as function of the waiting time . It 
can be seen that even if they could be fitted with power 
laws (with an exponent a ^ 0.38 ± 0.1 different from the 
exponent found for the Lenard- Jones case), the collapse 
of the different curves is not good. 

It is worth to make a comparison of our numerical re- 
sults with recent data measured by Pham et.al '27], for a 
system of colloids with depletion interaction. In their ex- 
periment, they measure the density-density correlation 
function in the aging regime for two system at differ- 
ent polymer concentration, corresponding respectively to 
conditions where the repulsive and the attractive glass 
are expected. To better compare with Pham et. al data, 
we calculated the density dynamical structure factor (ir- 
respective of the particle type) for the same g-vector of 
the cited experiment. Results are shown in Fig. |51 and 
Fig. IIQI where the T/ = 1.2 and Tf — 0.3 quenches are re- 
spectively plotted together with the experimental results 
from Ref. [23| (lower panel). Slight differences in the 
plateau values can be ascribed to differences in the ac- 
tual composition of the simulated (binary mixture) and 
experimental (slightly polydisperse system). The agree- 
ment between the numerical and the experimental data 
is impressive. 

Another important two-time observable is the mean 
squared displacement (MSD). In aging, it is defined as 
{\r{t — tyj) — r(tuj)P), where r(t) is the position of the 
particle at time t and the average (•) is performed over all 
the particles and the independent configurations. Here 
we focus on the total MSD, i.e. evaluated with no dis- 
tinction between particles of the two species. 

Fig. [nil shows the MSD for T/ = 1.2. In analogy 
with the case of the density correlators, at short times 
there is no significant iu,-dependence. At larger times 
the diffusive process is visible. When the waiting time 
becomes large, a plateau starts to emerge at intermediate 
times, between the ballistic (short time) and the diff'usive 
(long time) regime. The duration of the plateau grows 
as tw increases. As above, this behavior is similar to the 
one observed in previous simulations of a binary mixture 
of Lenard- Jones particles [2^ . 

This three time-region behavior (ballistic, plateau, dif- 
fusive) is generally explained in term of the so called cage 
effect. When repulsive interactions start to be dominant 
the system slows down since each particle starts to get 
trapped in a shell of first neighbor particles. The length 
of the plateau indicates the typical time that a parti- 
cle needs to escape this cage. The value of the plateau 
stands at the typical size of a cage, i.e. around 10% of 
the diameter in agreement with Lindemann's melting cri- 
terion [33 ■ For Tf = 0.3, the MSD shows a completely 
different behaviour, fFiglllb'l. Its evolution can again 
be divided in three time-regions, but now the short time 
(ballistic) and the long time (diffusive) regions bracket 
a long intermediate region where the MSD appears to 
grow as a power-law in time, with an effective exponent 
which decreases with i^,. Differently from the previous 
case, no evident sign of a plateau is present, except for a 



hint of a plateau at MSD = 0.001 which starts to develop 
for large waiting times for t — t^, « 0.5. This value is 
consistent with a displacement of the order of the width 
of the potential well, and can be interpreted as a result of 
particles bonding in the attractive well, i.e. as some sort 
of attractive cages. However, activated bond-breaking 
processes appear to destroy such a confinement on time 
scales longer than f — f„ w 1. 



IV. CONCLUSIONS 

In systems with short-ranged attractive potentials, an 
efficient competition between attraction and excluded 
volume generates a highly non trivial equilibrium dynam- 
ics and two kinetically distinct glasses. In the "attrac- 
tive glass" , the particle localization is controlled by the 
bonding distance, while in the "repulsive glass" localiza- 
tion is controlled by the neighbor location. This offers 
the unique possibility of studying two different aging dy- 
namics within the same model. In this article, we have 
realized this numerically for a model system which shows 
both types of disordered arrested structures, the square 
well potential. We have found that the aging dynamics is 
clearly different for the two glasses, both in the waiting 
time evolution of the one-time quantities, and in the time 
evolution of the two-times quantities. 

The results reported in this article refer to a simple 
model system which capture the essence of the physics 
introduced by the short-range attraction. The square 
well potential is indeed a satisfactory representation of 
the coarse-grained particle-particle potential in colloidal 
system, for example in the presence of depletion interac- 
tions. Of course, dynamical features of the solvent con- 
stituents and of the solvent-particle interactions may sig- 
nificantly affect the particle dynamics. In this respect, 
the experimental validation of the equilibrium results is 
very valuable. 

In particular, we have shown that the partial static 
structure factor evolves with i^, , showing opposite trends 
in the two cases, i.e. the first peak grows with t^, for the 
repulsive glass while decreases for the attractive one. We 
have also shown that the decay of correlation functions 
becomes slower and slower on increasing tw in both cases, 
but with noteworthy differences in the strength of the re- 
laxation (much more intense for the attractive glass case) 
and in the shape of the relaxation itself. The analysis of 
the MSD brings clear evidence that, while in the repul- 
sive case there is the emergence of a plateau, whose time 
duration increases significantly with t^^, , in the attractive 
glass case this does not happen. Indeed, in short-ranged 
attractive colloids, activated processes can be associated 
with thermal fluctuations of order uq, which allow parti- 
cles to escape from the bonds. These processes generate 
a finite bond lifetime and, at the same time, may desta- 
bilize the attractive glass. It is important to find out 
how the difference in the aging dynamics arising from 
the different localization mechanisms are affected by the 
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presence of activated processes. 

Finally, we have shown that a good agreement is found 
between experimental and numerical data, also in the 
aging dynamic. 
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FIG. 1: Square Well iso-diflusivity curves (i.e. the locus 
of points in the <j>, T plane where the normalized diffusion 
coefficient is constant) for three different values of D/Dq 
(-Do = ai,^T/m data from Ref.j2&])- The extrapolated ideal 
glass line is also reported (data from Ref.j^^). Equilibrium 
starting configurations at = 0.608 and Ti = 0.6 (filled cir- 
cle) are instantaneously quenched at f = to T/ = 0.3 and 
T/ = 1.2 (plus). 
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FIG. 3: First peak of the static structure factor at different 
waiting time t^j, 10^ (circles), 10^ (squares), 10^ (diamonds), 
10^ (triangles up) and 10* (triangles left). The quench is at 
Tf = 1.2. For long there is slow growth of the peaks of the 
structure factors. In the inset structure factors for the same 
tjn are displayed. 
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FIG. 4: Same as Fig. for T/ = 0.3. 



FIG. 2: Average potential energy U after the quench as func- 
tion of the time elapsed from the quench. The symbols are 
the simulation data and the continuous line are the fits. Data 
for t-u) < 1 are not shown since the kinetic energy has mot 
yet equilibrated to the final value. For Tf = 1.2 the fit is 
U/N = -2.77 -I- 0.01 In (t„) while for Tf = 0.3 we obtain 
[7/iV = -4.05 + 0.43 
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FIG. 5: (a) Density-density correlation function 

4>AA{q,tw,t — tw) after a quench at T/ = 1.2 for qas ~ 8.72. 
From right to left, longer and longer relaxation time, — 10^ 
(continuous line), 10^ (dotted line), 10'^ (short dashed line), 
10* (long dashed line), 10^ (dot dashed line), 
(b): same results but for gas — 6.37 
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FIG. 6: Same as Fig. but for the quench at Tf = 0.3. In 
the inset in (a), an enlargement of the plateau area is shown 
(see text for details). 
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FIG. 7: Scaling of the density-density correlation functions 
for Tf — 0.3. The waiting times are tw = 10^ (continuous 
line), 10^ (dotted line), 10^ (short dashed line), lO" (long 
dashed line), 10^ (dot dashed line). 

(a) best scaling for the region where 0.7 < (j)q < 0.9. 

(b) best scaling for the region where 0.35 < 4>q < 0.5. 
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FIG. 8: Scaling coefficient used to obtain the scaling in 
FiglTJi-b. The circle are the scaling coefficient for the case 
in Fig|7^, the square for|7j3. 
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(a) Simulations 
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FIG. 10: Upper panel: Same as in FigElfor Tf = 0.3. 
Lower panel: Again experimental data are taken from 
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FIG. 9: Upper panel: total density-density correlation func- 
tions for qRA = 2.93 at different waiting times for Tf — 1.2. 
From right to left, longer and longer relaxation time, — 
10\ 10^ 10^ 10", 10^ 

Lower Panel: Experimental data reproduced from [2^ . The 
unit of time are in seconds but time has been rescaled by the 
relative viscosity of the solvent and the wave vector. 



FIG. 11: (a) Mean square displacement (MSD) for the quench 
at Tf — 1.2 for increasing waiting times tw = 10^ (continuous 
line), 10^ (dotted line), 10^ (short dashed line), 10* (long 
dashed line), 10^ (dot dashed line) 
(b) As in (a) for Tf = 0.3 



